#!/usr/local/bin/Rscript

#
#  Barplot and pie chart from the output of annovar (.variant_function)
#
# Input looks like this (NB: No header, only first column actually used)
# -----------------------------------------------------------------------------
# > less *.variant_function
# upstream  FAM102B chr1    108904064       108904065       0       A       Merged-chr1-108904064-3
# UTR5      ATP5F1  chr1    111793500       111793501       0       A       Merged-chr1-111793500-5
# intronic  MAN1A2  chr1    117741372       117741373       0       A       Merged-chr1-117741372-2
# UTR5      WARS2   chr1    119484812       119484813       0       A       Merged-chr1-119484812-7

library(RColorBrewer)

anntable<- read.table('macs/mergepeaks/macs.mergeBed.annovar.variant_function', sep= '\t')
anncounts<- table(anntable$V1)
anncounts<- anncounts[order(anncounts)]
perc<- anncounts/sum(anncounts) * 100

options(warn= -1)
w= 8
h= w * 1
## bitmap('macs/mergepeaks/barplot_annotation.png', units= 'cm', width= w, height= h, res= 512, pointsize= 6.5)
png('macs/barplot_annotation.png')
par(mar= c(4, 10, 2, 4), las= 1, xpd= TRUE)
b<- barplot(anncounts, names.arg= names(anncounts), col= 'salmon', horiz= TRUE, main= 'Peak annotation by region', xlab= 'Number of peaks')
text(labels= paste(formatC(perc, digits= 2), '%', sep= ''), x= anncounts + max(anncounts)*0.05, y= b, adj= 0)
dev.off()
## system('scp macs/barplot_annotation.png berald01@143.65.172.155:~/Tritume/')

MAX_CTGS<- 4 ## Number of categories to include, everything else will be lumped into 'other' category 
if(length(anncounts) > MAX_CTGS){
    piecount<- c(rev(anncounts)[1:MAX_CTGS], sum(rev(anncounts)[(MAX_CTGS+1):length(anncounts)]))
    names(piecount)<- c(rev(names(anncounts))[1:MAX_CTGS], 'other')
    } else {
    piecount<- rev(anncounts)
    names(piecount)<- rev(anncounts)
}
w= 8
h= w * 1
## bitmap('macs/mergepeaks/pie_annotation.png', units= 'cm', width= w, height= h, res= 512, pointsize= 8)
png('macs/pie_annotation.png')
par(mar= c(2, 4, 1, 6), las= 1, xpd= TRUE)
pie(piecount, labels= paste(names(piecount), piecount), main= 'Peak annotation by region', col= brewer.pal(MAX_CTGS+1, "Dark2"))
dev.off()
## system('scp macs/pie_annotation.png berald01@143.65.172.155:~/Tritume/')